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Abstract 

Estimating the cosmological microwave background is of utmost importance for cosmology. However, its 
estimation from full-sky surveys such as WMAP or more recently Planck is challenging : CMB maps are generally 
estimated via the application of some source separation techniques which never prevent the final map from being 

■ contaminated with noise and foreground residuals. These spurious contaminations whether noise or foreground 
I residuals are well-known to be a plague for most cosmologically relevant tests or evaluations; this includes CMB 

^ ■ lensing reconstruction or non-Gaussian signatures search. Noise reduction is generally performed by applying a 

O I simple Wiener filter in spherical harmonics; however this does not account for the non-stationarity of the noise. 

■ Foreground contamination is usually tackled by masking the most intense residuals detected in the map, which 
makes CMB evaluation harder to perform. In this paper, we introduce a novel noise reduction framework coined 
LIW-Filtering for Linear Iterative Wavelet Filtering which is able to account for the noise spatial variability thanks 
to a wavelet-based modeling while keeping the highly desired linearity of the Wiener filter. We further show that 
the same filtering technique can effectively perform foreground contamination reduction thus providing a globally 
cleaner CMB map. Numerical results on simulated but realistic Planck data are provided. 
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p_i. I. Introduction 

Q . In mid 2009, ESA (European Spatial Agency) put in orbit the latest space observatory Planck to investigate the 

^ ' Cosmological Microwave Background (CMB). These data are of particular scientific importance as it will provide 
more insight into the understanding of the birth of our Universe. Most cosmological parameters can be derived 
from the study of these CMB data. After a series of successful CMB experiments (Archeops, Boomerang, Maxima, 

I _ COBE and WMAP ||T1), the Planck ESA mission is providing more accurate data which will be the reference full- 

J> ' sky high resolution CMB data for the next decades. More precisely, recovering useful scientific information requires 

. disentangling the CMB from the contribution of several astrophysical components namely Galactic emissions from 

^ I dust and synchrotron, Sunyaev-Zel'dovich (SZ) clusters, Free-Free emissions, CO emission to only name a few, see 

ff^ ■ 0. Classically, the CMB map is estimated via the application of some source separation methods. However, the 

^ ■ application of very sophisticated source separation methods (see Q, H, ||5l, ISl) does not prevent the final estimated 

y—^ map from being contaminated with various spurious components : i) instrumental noise and ii) foreground residuals. 

^ Dealing with these various sources of contaminations is of paramount importance for most cosmologically relevant 

• • tests or evaluations : this includes CMB lensing reconstruction I7| or non-Gaussian signatures search. Tackling the 

. ^ problem of noise and foreground residual reduction is therefore important and is at the heart of this paper. 

^: 

. a) Instrumental noise: The very specific scanning pattern in full-sky CMB experiments leads to an instrumental 

■ ■ ■ noise closely Gaussian that exhibits significant spatial variation : Planck noise is far from being stationary. As an 
illustration. Figure [T] displays the root mean square (RMS) map of Planck expected instrumental noise at 217GHz. 
Furthermore, while noise is insignificant at large scale, it is largely predominant at small scale (e.g. typically for 
i > 1500 for Planck). 

Most source separation methods applied so far to the Planck data are linear methods : the estimated CMB map can 
be expressed as a linear combination of all the pixels of the input observations. As such, the propagated instrumental 
noise in the CMB maps is also typically assumed non-stationary Gaussian. 

Noise can be a significant limitation for most cosmological tests or reconstructions. High noise level can hamper 
non-gaussianity tests as the noisy map is likely to be "more Gaussian" than the noise-free map. The classical 
approach in the field of cosmology generally consists in reducing the noise of the CMB via a Wiener filter in 
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spherical harmonics. If Ci stands for the theoretical power spectrum of the CMB and C" the power spectrum of 
the noise, the Wiener filter is defined as follows : 

V£,m>0; aL = ^-|^aL (D 

where {a£m}e,m {resp. {a^^}i^rn) stands for the spherical harmonics coefficients of the filtered map x (resp. the input 
map y). While this approach makes profit of some knowledge of the energy distribution of the noise in frequency, 
it does not account for its spatial variations or non-stationarity. It is very likely that such a non-stationarity may 
create global non-gaussian signatures at high frequency when such variations are sharp. Strong non-stationarities 
are known to have a strong impact on CMB lensing reconstruction (see [SJ), creating an undesired reconstruction 
bias. 

Another less traditional but commonly encountered noise reducing technique is the local Wiener filtering. It boils 
down to apply the Wiener filter in the Fourier space on patches in the pixel domain. The main drawback of this 
method is that while large patches should be favored to capture the non-stationarity of noise, they are enable to 
capture correlations at scales larger than the patch size. 

b) Foreground residual: The first evaluation of already sophisticated CMB -dedicated source separation meth- 
ods has been performed in [3|. Since, novel very effective methods have been proposed Needlet ILC [4|, Generalized 
ILC 19], SMICA ifTOll or GMCA [6|. However, due to the very high complexity of the data, none of these very 
effective methods is able to provide a CMB map that is guaranteed to be free of foreground contaminations. If 
the residual contamination is generally low at large scales ( e.g. for i < 500 for Planck), it is generally significant 
at higher frequencies. It is likely that the remaining contamination mainly originates from point sources or strong 
features in the galactic plane. Obviously the presence of these residuals highly biases any cosmological tests to be 
performed on contaminated CMB map. The usual approach consists in masking the known spurious pixels of the 
CMB prior to any post-processing (e.g. 20 to 40% of the sky are usually masked). However masking raises several 
issues : i) it limits the number of samples to which any test can be applied, limiting its statistical relevance and 
ii) masking generally makes any post-processing much more complicated to perform properly (e.g. inpainting of 
the mask is generally needed prior to any CMB lensing reconstruction ||3). Avoiding or at least limiting the extent 
of the mask to be applied should greatly help improving any cosmological test to be made on the estimated CMB 
map. Furthermore, masking does not prevent from the presence of residuals of undetected point sources which are 
likely be be numerous all over the sky. Further reducing the amount of foreground residuals should clearly be helpful. 

c) Contribution of this paper: The contribution of this paper is twofold : 

• It introduces a novel noise reduction framework that, in opposition to the classical spherical-harmonics based 
Wiener filter, accounts for the non-stationarity of the noise. This new method preserves the linearity of the 
filtering. In fact, linearity is crucial in this field as it makes the study of error propagation via Monte-Carlo 
simulations much more convenient to carry out. In this framework, we make use of a simple, but efficient, 
wavelet-based modeling of the noise that allows for the modeling of non-stationarities at different scales. This 
will be discussed in Section |lll] 

• We discuss an extension of this method to further estimate an approximate contribution of the foreground 
residuals per wavelet scales by means of RMS maps. Thereby contamination reduction could be tackled 
within the same filtering framework. 

The noise/contamination reduction problem, recast as a quadratically regularized least square problem, is then 
effectively solved using recently introduced algorithms based on proximal calculus. Extensive numerical experiments 
are provided that show the effectiveness of the proposed method in reducing noise in Section III-CI and foreground 
residuals in Section IIII-CI 

II. CMB MAP FILTERING WITH NON-STATIONARY NOISE 

A. CMB map filtering 

In a very general context and more precisely in the context of Planck, it is customary to assume that the input 
noise data, denotes by y in the following, is modeled as the linear combination of the sought after noise-free signal 
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Figure 1. Instrumental noise in Planck - simulated noise root mean square map at 217GHz in mK (Antenna temperature). 
X and the noise term n : 

y = X + n 

Decoirelation between x and n is also classically assumed. In what follows, we will assume that each of these 
signals is sampled on the sphere according to the Healpix ifTTI sampling grid. The sampled CMB signal x is 
assumed to be composed of N samples {x[k]}k=i,--- ,n- Recovering an estimate x of x can then be formulated as 
an inverse problem. According to the Bayesian way of solving inverse problems, some a priori knowledge about the 
signal of interest x has to be expressed. This inference framework is particularly well suited to the CMB denoising 
issue as strong prior assumptions can be formulated. Indeed, since |12!], it is well established that the cosmological 
microwave background can be approximated as an isotropic Gaussian random field with a high accuracy. This 
particularly means that it is fully characterized by its power spectrum Ci. 
More precisely, the CMB x can be defined by its expansion in spherical harmonics : 

oo i 

a; = ^ ^ aimYim (2) 

i=0 m=-i 

Assuming that the CMB is fully Gaussian with power spectrum Ce, its spherical harmonics {aemjern are indepen- 
dently distributed following a Gaussian distribution with mean and variance C£ : 

a^™~AA(0,Q) (3) 

The power spectrum of the CMB is generally estimated with the pseudo power spectrum defined as follows : 

Ci = 2£a_i a}^aem (4) 
m=-e 

where stands for the conjugate of a£m- The a priori probability distribution V of the CMB map x is a normal 
distribution with covariance matrix F^CF where F stands for the spherical harmonics basis and F^ its adjoint. 
The matrix C is diagonal with C = diag {Ce). This yields the following CMB prior probability distribution : 

V{x) oc exp -x^F^C-^Fx (5) 

In this section, we assume that the instrumental noise can be modeled as a non-stationary but uncorrected Gaussian 
field in the pixel domain such that : 

VA: = I,--- ,iV; n[A;] ~ (O, a|) 
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Written differently, n is a realization of a multivalued Gaussian variable with mean and covariance matrix 
S = diag (cjf • • • (T^) . The case of correlated noise will be discussed later on in this section. From the Gaussian 
modeling of n, the likelihood function is trivially obtained as follows : 

£(7/|a:,S) ocexp-(y-x)^S~^(y-x) (6) 

Following Bayes rule, the posterior distribution for x is the following : 

P{x\y) oc C{y\x,'S)V{x) (7) 



oc exp 



x'^T^C-^Tx + {y- xf S"^ {y 



(8) 



The CMB map can be estimated using the maximum a posteriori (MAP) estimator with the following analytical 
solution : 

X = Argmax^. (9) 
= {'^-^ +T^C-^Ty^'E-^y (10) 

The reader would here recognize a classical Wiener filter. Contrary to the classically used Wiener filter in spherical 
harmonics, the solution of Equation Q explicitly account for the non-stationarity of the noise. However, the 
computational evaluation of this Wiener-like solution raises some numerical issues : the noise covariance matrix 
is diagonal in the pixel domain while the CMB covariance matrix is diagonal in the spherical harmonics domain. 
The solution filter that appears in Equation ^ is formulated in the pixel domain : 

In this expression, the matrix T^C"^T is the inverse covariance matrix of the CMB in the pixel domain. This 
matrix is clearly not diagonal. Recalling that N"^ is of the order of lel3 for WMAP and 2el5 for Planck, storing 
and handling this matrix is way too unrealistic. As an alternative, the most effective approach in this situation 
would amount to solve the linear system of equations (S^^ + J^^C^^J^) x = using a conjugate gradient 

(CG - 1 13 1) for instance. However, in the next section, we will rather make use of recently introduced proximal 
algorithms. This optimization framework leads to algorithms which exhibits a lower speed of convergence than CG 
but with the gain of more flexibility to handle different filtering techniques or better constrain the solution. 

1 ) Iterative filtering: As said clearly in the previous paragraph, computing the MAP solution requires inverting 
a system of equations which turns to be intractable due to the large scale of Planck data. The most straightforward 
numerical solver is the well-known conjugate gradient. However, while this solver is known to be a fast and accurate 
numerical method for solving linear systems of equations, it lacks the flexibility of more modem optimization 
techniques. For this reason we rather opt for the more flexible optimization framework of proximal calculus (see 
fl4\ and references therein). 

a) A fiexible optimization framework - forward backward splitting: The problem in Equation (1101 ) can be 
written as minimization of the sum of two operators /i and /2 : 

X = Argmin^/i(x) + /2(x) (11) 
= x'^F^C-^Fx+iy-xf^^'^iy-x) (12) 

The function fi{x) = x^ F^Q>~^Fx is convex and lower semi-continuous and differentiable; the function /2(a:) = 
(y — x)^ (y — x) is strictly convex, continuous and differentiable. Such problem as therefore a unique solution 
(see mi ). 

If V/2(x) denotes the gradient of /2, we will also assumed that this gradient is L-Lipschitz. In the specific setting 
of interest in this article, V/2(x) = —2Tr^ {y — x) and L = 2||S~^||2. In this context, an effective algorithmic 
framework is the forward-backward splitting (EPS) technique defined in |fT4l|. In the paradigm of proximal calculus, 
any convex lower semi-continuous function (see |[T4l ) admits a well defined proximal operator : 

prox^^(z) = Argmin„7/('u) + ^\\z- u\\j^ (13) 
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In the FBS framework, the solution to the problem in Equation (fTTI ) is computed via an iterative fixed point 
algorithm such that at iteration k : 

= prox^^^ (x'^ - tV/sIx'^)) (14) 

= prox^^^ (x'= + 27S-i (y-x'^)) (15) 

As proved in f\A\, the convergence of the above fixed point algorithm is guaranteed under the condition ^ < j-; ; 
note that in the setting of CMB map recovery with uncorrelated non-stationary noise, L = 1/ min^ a^. 

b) The proximal operator of fi: In the problem of interest in this paper, /i(x) = C~^Fx. As 

this function is continuous and differentiable, it is a fortiori lower semi-continuous. Following the definition in 
Equation ([T3]) . the proximal operator of /i is defined as follows : 

prox^^^ (z) = Argmin„7u^J"^C"^J'n + ^ll-z - u\\\ (16) 
The solution to this problem has a simple closed-form expression : 

prox^x-.F«c-.Fx = (/ + 27^^C-ij-)"'z (17) 
= + Fz (18) 

The proximal operator of /i is a mere Wiener filter in the spherical harmonics space. 

To sum everything up, the resulting iterative Wiener-like filtering is as follows : 



1. Initialization : Set the number of iterations /max, tine step size 7 = miiifc ol and the starting point x 

2. Main loop : apply iteratively the following steps at each iteration k : 

a - Forward step of the FPS : 

u^x'' + 27S-^ (y - x'') 
b - Backward step of the FPS : 

3. Stop when k = /max- 



c) Computational cost: The computational cost of this iterative Wiener algorithm is particularly low. The first 
step (2 - a) only requires entry wise multiplications and additions of vectors. It is therefore of the order of 0{N). 
The second step (2 - b) is by far the most computationally demanding; it requires entrywise multiplications and 
additions of vectors of the order of 0(£max) (where generally ^max — 3000) and single forward/backward spherical 
harmonics transforms. In the numerical experiments of Section ITl-CI and ITlI-CI the starting point x^^'^ will be chosen 
as having entries equal to zero. The total number of iterations will be no greater than /max = 250. 

d) Linearity of the iterative Wiener filter: The proposed iterative filtering technique preserves the linearity of 
the overall CMB map processing. In fact, in every step of the proposed algorithm, each output depends linearly 
on their inputs. This property is essential as it allows for a simple way to study the propagation of errors via 
Monte-Carlo simulations; each simulation then has to undergo exactly the same processing step described above. 

B. The case of correlated noise 

In most CMB surveys, such as WMAP or Planck to only name two, the instruments have been thoroughly 
calibrated on ground thus meaning that the noise level of these instruments is accurately known. Assuming a 
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perfect calibration of the data, the non-stationarities of instrumental noise mainly originate from the non-uniform 
scanning of the sky. The lower the number of scanning time in one area is, the higher the noise level is and 
vice versa. Figure \T\ shows the RMS map of noise at channel 217GHz; mathematically this is no more than the 
square root of the diagonal of S. However this is also well known that, contrary to what has been assumed in the 
previous section, the noise is far from being uncorrelated. There are mainly two origins for the correlation of noise : 

• The way the sky is scanned by Planck makes the noise correlated along the scanning direction. This means 
that the noise is correlated along elongated patterns in the pixel domain; for numerical reasons, this can hardly 
be accounted for in the pixelwise covariance matrix S. 

• Most modern CMB-dedicated source separation methods Needlet ILC |4], GeneraUzed ILC [9], SMICA IlOl 
or GMCA |6] rely on a local analysis of the data in the wavelet or needlet domain. Furthermore, the final 
CMB map is generally obtained as the mixture of observations that do not share the same resolution. This 
means that the noise that finally contaminates the CMB is, in the pixel domain, correlated. 

It then seems natural to model non-stationary correlated noise : i) in the wavelet domain to capture its correlation 
and ii) locally to capture its non-stationarity. 

1) Wavelets - the basics: Wavelets are well known to be the tool of choice to analyze non-stationary signals 
|[T5l . When considering spherical data, there exist several implementations of wavelets on the sphere. In this paper, 
we will make use of the isotropic undecimated wavelet transform (lUWT) introduced in |[T6l . The lUWT can be 
seen as an extension of the celebrated a trous algorithm to spherical data (see the seminal book |[T5l and references 
therein). 

As any wavelet transform, the lUWT is first defined by a scaling function 4'i^{0, ip) defined by a fixed frequency 
cut-off ic- As emphasized in |[T6l . this scaling function is built so that it is azimuth-invariant on the sphere {i.e. it 
does not depend on ip). Classically, each rescaled version of the scaling function (j)ej^ is associated to a low-pass filter 

2-) 

hj or equivalently Hj{i) in spherical harmonics. Traditionally, the wavelet-based multiscale analysis is performed 
by computing successive smooth approximations of x. These smooth approximations are obtained by convolving 
X with dyadically rescaled versions of the scaling function : 

Vj = 0, • • • , J; Cj = hj -k X 

where J is the total number of scales. The so-called wavelet function at scale j, ipea{d,^)^ is defined as the 
difference of two consecutive rescaled scaling functions : 

V-ic {9, p) = 0_«c {9, (f) - (jjia {9, ip) 

23 23^1 23 

The wavelet function at scale j can be equivalently associated to a high-pass filter gj or equivalently Gj{£) = 
1 — Hj{£) in spherical harmonics. The wavelet coefficients at scale j are then computed by convolving the smooth 
approximation of x at scale j with the j-th wavelet function '>ph^{9,p). 

23 

Vj = 0, - • • , J; Wj+i = gj-kCj 

In spherical harmonics, this amounts to a simple multiplication of the signal's spherical harmonic coefficients with 
the wavelet filter Gj. 

The decomposition being redundant, several strategies can be designed to reconstruct x from the wavelet coefficients 
and the coarse resolution cj. The simplest is the one advocated by the a trous algorithm : 

J 

x = cj + y^^wj 

In the following, noise modeling will be performed in each wavelet scale; i.e. on each set of coefficients {wj}j=i^... j. 
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2 ) Wavelet-based statistical modeling of correlated noise: As emphasized previously, one elegant way to model 
the instrumental noise is to consider the distribution of its expansion coefficients in the spherical wavelet domain. 
The basic idea consists in assuming that the wavelet coefficients of noise are approximately decorrelated. This idea 
takes its roots in the field of multiscale statistical modeling : it has been shown that wavelets exhibit (almost)- 
decorrelating properties for some classes of non-stationary stochastic processes. To give only one example, this 
is the case for fractional Brownian motion |17| which has been extensively used to model power-law following 
stochastic processes. In the following, w^[k] will denote the wavelet coefficient of the noise term n in scale j at 
pixel k with the convention : j = 1 corresponds to the finest scale and j = J — 1 to the coarse resolution. From 
this definition, we will assume that these wavelet coefficients verify : 

E{w][k]w][i]}=a;;[kf6k,i (19) 

where cr"[fc]^ stands for the local variance of the noise at pixel k and 5^ j for the kronecker delta function. This 
entails that the covariance matrix of the noise in wavelet scale j is diagonal and equal to : 

S„j=diag(cT,"[fc]2) (20) 

The local standard deviation of the noise, 5]„ j, has to be estimated from the data beforehand. In the numerical 
experiments in Section III-C[ these parameters are estimated from a single noise realization. Within Planck, it is 
common to compute the so-called jack-knife map defined as the difference between consecutive scanning rings; 
provided that the instrument calibration do not vary from one ring to the other, it gives a noise realization. 
Assuming that the local variances cj"[fc]^ of the noise vary slowly from one pixel to another, it can be approximated 
at pixel k by the empirical estimate of the variance in a fixed surrounding neighborhood; in the following this 
neighborhood will be a patch of size \/P x \/P centered about the pixel k : 

[fc]' = p E (21) 

Where J\f[k] stands for the indices of the pixels that compose the patch that surrounds the pixel k. Note that the 
patch size ^/P will highly depend on how fast the noise variance varies spatially. This particularly means that 
"v/P should be smaller in the finest scale and larger when j increases. Following the natural dyadic scaling of the 
wavelet transform, the patch size in scale j will scale as follows : Pj = Pi where Pi stands for the patch size 
at the finest scale. 



C. Numerical experiments 

In this section, we particularly focus on the noise reduction aspect of the proposed method. The reference 
denoising method in the field of astrophysics and more specifically CMB data analysis is the global Wiener filter 
applied in spherical harmonics. In opposition to this classical filtering technique, the proposed iterative filtering 
approach accounts for the non-stationarity of the noise; thereby, it should particularly provide better solution that 
global Wiener. 

As detailed in the previous section, the modeling of the noise contribution is performed in the wavelet domain : 
the noise is assumed to be non-stationary but decorrelated at each wavelet scale j. In each wavelet scale, the noise 
covariance matrix is estimated locally from a single noise realization. The noise covariance per wavelet scale is 
evaluated locally on patches of size A^^^'^^ where, by convention, j = 1 is defined as the finest wavelet scale. As 
an illustration. Figure [8] displays the estimated noise RMS map at scale j = 1. 

The global Wiener solution is computed by filtering the original CMB map y with the linear filter defined in 
spherical harmonics by Equation The CMB power spectrum is a WMAP7 best-fit estimate; the noise power 
spectrum is estimated via the pseudo power spectrum of the noise realization. 

Iterative Wiener has been applied to the original noisy CMB map y. The single parameter to be tuned is the 
number of iterations /max; in the following experiments, it has been fixed to 250. Figure |2] shows the noisy map, 
the LIW-Filtering solution and the noise realization used to estimate the noise local variance in the finest wavelet 
scale. Figure [3] displays the residual maps of the global Wiener and LIW-Filtering solution. As expected the non- 
stationarity footprint of the noise appears clearly on the residual LIW-Filtering map where more noise seems to 
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have been removed. 

Figure |4] features the power spectra of the original and filtered maps. We would like to point out that these power 
spectra have been computed from the raw maps; no masking has been applied. This particularly means that these 
estimates are likely to be biased by strong foreground from the galactic plane. The first observation is that the 
global Wiener filter does not fully remove the noise at high i contrary to the iterative filtering technique. This 
can be explained by the presence of significant non-stationarities at high I which can produce a substantial bias at 
high frequency. Interestingly, the iterative technique succeed in largely reducing noise at high i thus supporting the 
significant contribution of the noise non-stationarities in producing a bias on the power spectrum. More technically, 
it is worth noting that the power spectrum of a stochastic process is properly defined as long as this process is wide 
sense stationary |[T8l . Obviously, the non-stationarities of the noise make it depart from this assumption; thereby 
the power spectrum has no rigorous statistical meaning in this context. Said differently, the power spectrum of the 
noise C" does not contain all the information that a non-stationary noise does in the pixel domain. 
It is well known that the Wiener filter provides a CMB map that has a biased power spectrum. Thereby the Wiener 
filter - in red in Figure |4] - yield a power spectrum that is negatively biased at high I when the noise highly 
predominates upon the CMB - typically for i > 2500. 

Figure |5] shows the power spectrum of the residual : x — x. This quantity is particularly well suited to visualize 
the power discrepancy in frequency between the original CMB map x and its estimate x. As expected, the global 
Wiener filtering reduces the contribution of noise at high i (see the blue line in Figure |5]). Still, the filtered map 
largely departs from the original map. Interestingly, the proposed iterative filtering technique (in red in Figure [5]) 
largely outperforms the classical Wiener filter from i = 900 to ^ = 3000. More precisely, the residual power is 
reduced by up to an order of magnitude for £ > 2500. 

A more complete measure of discrepancy between the original map x and one of its estimate x consists in evaluating 
the mean squared error (MSE) in the wavelet domain. To that end. Figure [6] displays the MSE of different estimators 
of the CMB map in 5 wavelet scales in logarithmic scale. The normalized MSE between some CMB map estimate 
X and the "true" CMB map is computed as follows : .rt.i2 Again, this figure shows that the iterative 

technique which directly account for the non-stationarity of the noise, performs better than the global Wiener filter. 
The impact of the noise tends to vanish for larger scales {typically for i < 1500). 

To better illustrate the differences between the global Wiener filter and its iterative counterparts. Figure [7] displays 
the MSE of the different estimators of the CMB map which has been computed in the wavelet domain per bands of 
latitude of width 10°. By convention, the galactic plane is centered around the latitude 0°. The figure on the right 
displays the MSE which has been evaluated from the finest wavelet scale j = 1; this corresponds to an analysis 
window that ranges from i = 1500 to ^ = 3000. Again, taken globally, one can point out the clear improvement 
between the global Wiener filter and its iterative version. As expected the Wiener filter is optimal in the sense of 
the MSE. More interestingly, the MSE of the two filtering techniques seems to be rather constant for high latitudes 
and sharply increases in the range [—15°, 15°]. As already emphasized, the impact of any of these linear filtering 
techniques is the more significant when the signal to noise ratio (SNR) decrease. This is obviously the case at 
smaller scales. At larger scales, the SNR increases and the impact of the filtering is much less significant. This is 
exactly what can be observed on the plot on the right of Figure [T] : the MSE of any of the filtered CMB maps at 
scale j = 2, which grossly correspond to an analysis window ranging from £ = 750 to ^ = 1500, are approximately 
the same. 



As a brief conclusion for these experimental results, it appears clearly that the non-stationarity of the noise has 
an impact on the CMB estimate. The most classical noise reducing technique in the field, i.e. the global Wiener 
filter, only relies on the power spectra of the CMB and the noise C". Obviously, the behavior of the noise power 
spectrum contains almost no information about its spatial behavior. These experiments have shown that the noise 
non-stationary has a clear impact on the CMB. Further accounting for the spatial variations of the noise variance 
helps improving the noise reduction process by a dramatic amount. 

III. From noise reduction to contamination reduction 

As stated in the introduction in Section Jl one of the purposes of this paper is to tackle the problem of post- 
separation contamination reduction. In this setting, contamination generally means instrumental noise and residuals 
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Figure 2. CMB maps - Top-left : Noisy image. Top-right : LIW-Filtering solution. Bottom-Left : Residual map (true map - noisy map) 
in the finest wavelet scale. Bottom-Right : Residual map (true map - LIW-Filtering solution) in the finest wavelet scale. Units in mK. 




Figure 3. Removed noise map (i.e. noisy data - filtered data) in the finest wavelet scale - Left : from the global Wiener filter. Right : 
from the LIW-Filtered solution. These maps correspond to what has been filtered from the noisy map. The amount of noise removed 
by LIW-Filtering is more realistic. Units in mK. 



of galactic foregrounds. In the case of Planck, CMB maps are generally obtained by applying sources separation 
methods. Most of these methods provide a solution that can be formulated as a linear combination of the original 
Planck multichannel observations. Therefore, it is relevant to assume that these contaminants essentially contribute 
additively to the CMB map : 

y = x + f + n (22) 
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Figure 4. Power spectra. In abscissa : spherical harmonics coefficient. In ordinate : CMB power in mK^. Dashed line show the power 
spectra of the noise for the different maps. 
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Figure 5. Residual power spectra. In abscissa : spherical harmonics coefficient. In ordinate : CMB power in mK^. Dashed line show the 
power spectra of the noise for the different maps. 

where / stands for the foreground contamination and n for the instrumental noise. Furthermore, it is also natural 
to assume that these three components x, f and n are mutually uncorrected. The filtering technique introduced in 
the previous sections mainly relies on the fact that the contaminants are characterized by their covariance matrix 
S. The questions that we aim at tackling in this paragraph is : How far do we know XI and how can it be estimated 
from the data ? 
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MSE per wavelet scale 




Figure 6. MSE per wavelet scale. In abscissa : wavelet scale - corresponds to j — 1. In ordinate : normalized mean squared error. 
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Figure 7. MSE per latitude at wavelet scale j ~ 1 and j = 2. In abscissa : latitude. The galactic plane corresponds to 0°. In ordinate : 
normalized mean squared error. 



A. Foreground residuals modeling 

Unlike noise, the foreground residual are generally non-Gaussian, with the exception of the CIB (see |fT9l ). We 
further restrict the modeling of residuals to their second order statistics thus opting for a Gaussian approximation 
of their contribution. The contribution of foreground residuals are obviously not known in advance; this means 
that their covariance matrix has to be estimated from the data directly. From the uncorrelation of the different 
components, the additive contamination model in Equation (l22l) then reads at the level of second statistics : 



^ — + + Yin 

In this equation, the correlation of foreground pixels implies that 5]y is non-diagonal. Similarly to the noise, second 
order dependencies of foreground pixels are way too complex to be modeled in the pixel domain. 
So as to capture the correlation of foreground pixels, a natural and simple strategy boils down to adopting the 
wavelet-based statistical modeling used to model correlated noise in the previous Section III-B2I To this end, Wj [k] 
will denote the wavelet coefficient of / in scale j at pixel k with the convention : j = 1 corresponds to the finest 
scale and j = J — 1 to the coarse resolution. From this definition, we will assume that these wavelet coefficients 
verify : 

K[w^.[k]w^.[i]] =af^[kf6k,, (23) 

where stands for the local variance of the foreground at pixel k. The covariance matrix of the foreground 

residual in wavelet scale j is diagonal and equal to : 



=diag(cjJ[A;]^ 



(24) 



In practical situations, Hfj has to be estimated from the data y themselves. Similarly to the modeling of noise, 
the local noise variance of the foreground will be approximated by the empirical estimate of the variance in a 
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surrounding neighborhood of size \/P x \/P centered about the pixel k 



(25) 



B. Contamination filtering 

In the remaining of this paper, we will generally use the term contamination for the contribution of both the 
foreground residuals and the instrumental noise. In the following, as suggested previously, we also opt for exactly 
the same model to capture the pixel dependencies of the instrumental noise. Therefore, we will now define locally 
at pixel k and scale j the variance of contamination (noise and foreground) as (t^[A;]^ = (5"J[A;]^ + (t"[/c]^. The joint 
contribution of noise and foreground residuals will be evaluated jointly from the data. 

From the estimation of the estimation of the local variance a^-[k\^; one has to recall that the contaminants' 
variance is not directly accessible but is equal to : 



(26) 



where (t| stands for the variance of the CMB and [x]+ = max(0, x). Let use note that, from the stationarity of 
the CMB map, crj does not explicitly depend on k\ this holds as long as the patch size at scale j is large enough 
compared to the CMB fluctuations (i.e. correlation length) in the j-th wavelet scale. The value of the CMB variance 
at scale j-th is directly related to the CMB power within this scale; its value can be derived from the CMB power 
spectrum Ci and the wavelet analysis filters. To that end, let us denote by ^ the spherical harmonics filter from 
which the wavelet scale j is defined. As emphasized in Section |I] the spherical harmonics coefficients of the CMB 
map should be mutually uncorrected. This particularly entails that the power of the CMB in the j-th wavelet scale 
is exactly : 

It follows that the variance of the contamination at pixel k is approximately equal to : 



(27) 



Thereby, from the a priori knowledge of the CMB power spectrum Ci, an approximate contribution of the residual 
can be evaluated at each pixel and each wavelet scale. The reader has to keep in mind that the measurement of this 
contribution is only approximate and has to be considered under the light of the following assumptions : 1) the 
contribution of the residual is modeled as a non-stationary Gaussian field where the covariance matrix is diagonal 
in each wavelet scales and 2) its variance being evaluated on small patches of size \/P x \/P, it is assumed to be 
approximately stationary in a given analysis patches. However, even with these approximations, our model offers 
a much more realistic modeling of the complexity of the data than the traditional Wiener method. Improvements 
of this modeling will be discussed in Section IIV-AI 

Let "Sf j = dmg{aj[k]'^). Then, in the algorithm introduced in Section |lll the noise matrix S is substituted with 
the following covariance matrix : 




























w 



where W stands for the isotropic wavelet transform and its adjoint. Along with the iterative technique 
introduced in Section JIl the resulting wavelet-based filtering will be coined LIW-Filtering for Linear Iterative 
Wiener Filtering in the following. 
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on line processing : 




Figure 8. Estimated noise covariance map at wavelet scale j = 1 : this covariance map has been estimated from a single realization of 
the simulated Planck noise that contaminates the CMB estimated by GMCA. Units in mK^. 

C. Contamination reduction on Planck simulated data 

As emphasized in Section |llll the wavelet-based noise/contamination variance modeUng used so far makes it 
possible to account for the contribution of the foreground contamination in the filtering process. Estimating the 
contribution of the foreground contamination, as detailed in Section JIIJ makes the reduction of contamination 
possible. In the following experiments, the contamination, which includes the contribution of the noise and an 
approximate contribution of the foreground residuals, has been modeled in the wavelet domain 1161 using 6 scales. 
Within each wavelet scale, with the exception of the coarsest scale to which no filtering is applied, the contribution 
of the contamination is assumed to be Gaussian and locally stationary on patches of size {4-'^^}_^._^ ^. Figure |9] 
features the estimated variance map in logarithmic scale. As a remark, if the morphology of this contaminant 
variance map resembles at first glance the noise-only variance map, the dynamic range of the former is clearly 
larger by 3 orders of magnitudes. These differences are mostly significant on the galactic plane and on local features 
that are likely point sources or residuals of SZ clusters. 

The proposed Wiener based iterative method has been applied on the aforementioned Planck simulated data. 
Figure [To] shows the difference map between the noise-only LIW-Filtering solution and the proposed contamination- 
aware LIW-Filtering solution in the finest wavelet scale. It appears clearly that the difference between the map are 
mainly concentrated on the galactic plane and on likely point sources at high latitude. 

Figure [TT] displays the power spectra of the original and filtered CMB maps along with the theoretical power 
spectrum of the simulated map in black dash-dotted line. At first glance, accounting for the contribution of the 
foreground residuals, even approximately, yields improvements from I = 750. Perhaps more interesting. Figure [12] 
shows the power spectrum of the residual maps x — x: the contaminant-aware filtering technique clearly outperforms 
the previous methods from i = 250 to ^ = 3000 with a dramatic improvement in the range [750, 1500]. 
Figure [13] displays the wavelet-based MSE we described in the previous paragraph. The results of the contaminant- 
aware filtering technique have been further added to the plot in Figure [6] Let us recall that the noise reducing 
techniques described in the previous section mainly help reducing the noise in the noise-dominant regime. This 
particularly explained why the MSE remains almost unaltered in wavelet scales j > 1. Interestingly, accounting for 
an approximate contribution of the foregrounds helps reducing their impact in larger scales, typically for j < 3. 
Figure [14] presents the MSE computed at different latitudes in the finest wavelet scale j = 1 - plot on the left - 
and j = 2 -plot on the right. Accounting for the contaminant contribution improves the MSE of the filtered map at 
all latitudes. Interestingly, the MSE is decreased by a large amount - a factor of approximately 5 - on the galactic 
(latitude 0°). The same holds in a slighter amount in wavelet scale j = 2. 

Assuming that the Planck instrumental noise is Gaussian is widely considered as a reasonable assumption. This is 
obviously not the case for foreground residuals the presence of which is likely to largely distort the search for non- 
gaussian features in the CMB. Thereby, reducing the amount of contaminant should help preventing non-gaussianity 
test from the non-Gaussian impact of foreground residuals. In the field of CMB non-gaussianity evaluation, a 
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classical technique boils down to measuring higher order statistics in the wavelet domain |[20l . One of the statistics 
of choice is the statistics of order 4, denoted K4, a.k.a. the kurtosis. The kurtosis of the wavelet coefficients of 
the estimated CMB map has been evaluated in 5 scales; non-gaussianity test results are shown in Figure [15] The 
first observation is that the proposed method helps reducing non-gaussian contamination even when only the noise 
is modeled as a spurious component. This suggests that the non-stationary behavior of the noise can generate 
undesired non-gaussian signatures. Thereby, accounting for the non-stationary behavior of the noise in the noise 
reduction process decreases the kurtosis of the filtered map by an order of magnitude in the finest wavelet scale. 
The red dashed line in Figure [15] shows a clear reduction of the kurtosis of the filtered map when the contribution 
of the foreground residuals is modeled. More precisely, the value of K4 is grossly reduced by 4 orders of magnitude 
for j = 1, 2 orders of magnitude for j = 2 and 1 order of magnitude for j = 3. While noise-only filtering yields 
map that departs from the Gaussian assumption for J < 2, the contamination-aware filtered map is now compatible 
with the Gaussian assumption in all wavelet scales with K4 < 1. 

Figure [16] presents more detailed non-gaussianity tests; the three plots that this figure displays show the evolution of 
K4, in the three finest wavelet scales j = 1, • • • ,3, which has been evaluated in bands of latitude of width 10°. The 
galactic plane is centered about the point 0°. The non-Gaussian contamination mainly originates from the presence 
of galactic foreground residuals which are mainly concentrated on the galactic plane. This explains the very large 
value of in the range of latitude [—15°, 15°]. The second well-known source of non-Gaussian signatures at 
spatial high frequency are the radio and infrared point sources. These foregrounds are roughly uniformly scattered 
all over the sky. Point sources are likely to be the predominant source of non-Gaussian signatures at high latitude 
thus explaining the relatively large value of K4 in this range of latitudes. As shown in the top-left plot of Figure [T6l 
further filtering the contamination (see the red dashed line) yields a dramatic reduction of K4 up to 5 orders of 
magnitude on the galactic plane. The filtered CMB map is likely to be compatible with the Gaussian assumption at 
every latitude, even on the galactic plane. In the second wavelet scale j = 2, the same reduction of non-Gaussian 
signatures is also visible by a lesser extent, mainly on the galactic plane in the range [—15°, 15°]. In the third wavelet 
scale, the proposed filtering is likely to reduce the non-Gaussian signatures in the galactic plane; in this range of 
scales (i.e. the third wavelet scale roughly corresponds to an analysis window that spans the range [375, 750]), the 
CMB is largely predominant which entails that : i) the input map is already quite compatible with the Gaussian 
assumption and ii) the impact of the filtering may be less visible when non-Gaussianity is measured. 



on line processing 




Figure 9. Estimated contamination covariance map at wavelet scale j — 1 : this covariance map has been estimated from the CMB 
estimated by GMCA. Units in mK^. 



IV. Discussion, prospects and conclusion 

A. Discussion and prospects 

a ) Filtering the maps, for what use ?: Important cosmological tests and evaluations are performed on estimated 
CMB maps; this is particularly the case for CMB lensing reconstruction [7J and non-Gaussian signature detection 
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Figure 10. CMB and residual maps Difference map between the noise-only LIW-Filtering solution and the contamination-aware LIW- 
Filtering solution in the finest wavelet scale. Units in mK. 
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Figure 11. Power spectra.. In abscissa : spherical harmonics coefficient. In ordinate : CMB power in mK^. Dashed line show the power 
spectra of the noise for the different maps. 

to only name two. These tests are particularly sensitive to spurious components that contaminate the estimated 
CMB map, whether it is noise or foreground residuals. This is therefore crucial to estimate a CMB map that is the 
less contaminated possible. The second point we would like to emphasize is that the aforementioned cosmological 
tests or evaluations are generally performed on the full-sky data; prior to any post-processing the estimated map is 
masked to get rid of the impact of galactic foregrounds and point sources. This has two major consequences : i) it 
limits the size of the samples to which these tests/evaluations are performed thus increasing the statistical variance 
of these tests and ii) any of these tests has to dealing with this mask (e.g. via an inpainting procedure in |7 |) thus 
making the analysis of the CMB map generally more complex. 

As shown in the previous experiments, the proposed filtering clearly limits the impact of noise by taking into account 
its non-stationary behavior. Furthermore, results are presented that clearly show that modeling the contribution 
of foreground residuals make contamination reduction effective even on the galactic plane. This has two major 
consequences : i) it makes it possible the use of a much smaller mask prior to any analysis of the estimated CMB 
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Figure 12. Residual power spectra. In abscissa : spherical harmonics coefficient. In ordinate : CMB power in mK^. Dashed line show 
the power spectra of the noise for the different maps. 
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Figure 13. Normalized MSE per wavelet scale. In abscissa : wavelet scale - corresponds to j = 1. In ordinate : normalized mean 
squared error. 



map and ii) it helps reducing the non-Gaussian features that originate from the presence of foreground residuals. 

b) Beyond the proposed methods: Whether noise or foreground residual is at stake, the central assumption that 
is at the heart of the modeling of contamination is the decorrelation hypothesis we made in Section |llll It is assumed 
that the contaminants have potentially non-stationary but decorrelated samples in each wavelet scales. The basic 
idea that supports this assumption is the well-known decorrelating property of wavelets bases. If this assumption 
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Figure 14. Normalized MSE per latitude at wavelet scale j ~ 1 and j — 2. In abscissa : latitude. The galactic plane corresponds to 0° 
In ordinate : normalized mean squared error. 
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Figure 15. Kurtosis K4 per wavelet scale. In abscissa : wavelet scale - corresponds to j = 1. In ordinate : value of the kurtosis K4. 



almost holds for rattier tlieoretical stocliastic processes (e.g. fractional Brownian motion ifTTl to only name one), 
this is only roughly approximate for more general signals such as galactic foreground or complex correlated noise. 
The basic idea behind this assumption is that the multiscale analysis basis is able to capture the correlation 
morphology of the component to be modeled. In the context of Planck, it is well understood that instrumental noise 
will be highly correlated along the scanning direction. This means that the correlation of the Planck instrumental 
noise be better modeled as a decorrelated stochastic process in a multiscale signal representation that best represent 
elongated structures such as the Curvelet tight frame [16|. We can extrapolate the same argument to the modeling 
of foreground residuals. 

The contamination modeling used so far in this article makes profit of the decorrelation assumption; it particularly 
help simplifying the filtering process by only requiring the handling of diagonal matrices (i.e. root mean squared 
maps). Departing from the decorrelation assumption will largely increase the complexity of the proposed filtering 
techniques. However, a straightforward way of extending the proposed methods is to choose multiscale signal 
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Figure 16. Kurtosis K4 per latitude at wavelet scale j = 1 and j — 2. In abscissa : latitude. The galactic plane corresponds to 0°. In 
ordinate : value of the kurtosis K4. 

representations which are better adapted to the morphology or structure of the signal to be modeled. 

c) Potential impact of contamination reduction on CMB non-gaussianities: It is important to wonder whether 
the contamination reduction filtering introduced in Section Hill might affect or reduce potential CMB-related non- 
Gaussianities. One important point is that the way the local variance of the contamination is estimated in Section Hill 
leads, in practice, to an estimation of foreground contaminations that largely exceeds the average CMB power in 
each wavelet scale. Therefore, the filtering should have a much lesser impact on non-gaussianities whose magnitude 
is of the order or lower than the magnitude of the CMB. It is worth pointing out that the rule in Equation l27l to 
estimate the local contamination variance can also be chosen differently to only account for non-gaussianities 
that are guaranteed to largely exceed the magnitude of the CMB thus avoiding for sure any effect on potential 
CMB-related non-gaussianities. Future work will focus on evaluating the potential impact of the proposed filtering 
technique on CMB lensing reconstruction and non-gaussianity detection. 

B. Conclusion 

Cosmological microwave background maps estimated from full-sky surveys such as WMAP or more recently 
Planck generally suffers from various sources of contaminations : i) instrumental noise is generally non-stationary 
which may generate non-gaussian signatures and ii) foreground residuals generally remain even after the application 
of state-of-the-art source separation methods. In this context, the most classical denoising technique, aka. global 
Wiener filter; despite its simplicity, it is not able to account for the non-stationarity of the instrumental noise. 
To that end, we introduce a novel noise reduction technique coined LIW-Filtering for Linear Iterative Wavelet 
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Filtering which combines the Unearity of the global Wiener filter while accounting for the potential non-stationarity 

of the noise. In this framework, the noise is modeled as a non-stationary but decorrelated process in the wavelet 
domain. The denoising problem then takes the very simple form of a quadratically regularized least square where 
the noise covariance matrix is diagonal in the wavelet domain and the signal covariance matrix is also diagonal 
in the spherical harmonic domain (i.e. CMB power spectrum). The solution is computed using recently introduced 
proximal algorithms. When noise reduction is at stake, we showed that the proposed iterative technique succeeds in 
reducing the mean squared error (MSE) of the filtered solution with respect to the global Wiener filter classically 
applied in the field. Moreover, the modehng/estimation framework introduced in this articles makes the reduction of 
foreground contamination possible. Similarly to the non-stationary instrumental noise, foreground contamination are 
modeled as a non-stationary but decorrelated process in the wavelet domain. Numerical experiments show that : i) 
the MSE of the filtered map is improved; specifically on the galactic plane and 2) very interestingly, non-Gaussian 
signatures are also dramatically reduced. This particularly provides arguments supporting the application of the 
proposed filtering technique as post-processing step to be applied to the CMB with the crucial aim of reducing 
noise and perhaps more importantly foreground contamination. It is also important to point out that, similarly to 
the global Wiener filter, LIW-Filtering is - at this acronym indicates - a linear filtering technique. This means that 
LIW-Filtering is also relevant when studying errors and their propagation on the estimated CMB map via Monte- 
Carlo simulations are unavoidable. 

Future work will focus on refining the modeling of noise and foreground residuals. Without losing the simplicity 
of this approach, we will more specifically focus on studying more adapted signal representation to better model 
the noise/foreground contamination spatial behavior. 

The developed IDL code will be released with the next version of ISAP (Interactive Sparse astronomical data 
Analysis Packages) via the web site: 

http://jstarck;.free.fr/isap. html 
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